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I.  SUMMARY 

The  research  performed  under  the  contract  during  the  period 
1  November  1981  through  30  April  1982  can  be  divided  into  two 
main  topics;  coupling  of  surface  waves  in  laterally  inhomogeneous 
source  regions  to  teleseismic  propagation  paths,  and  a  review  of 
the  theory  and  application  of  synthetic  seismograms. 

In  Section  II,  the  dimensions  of  the  cylindrical  source 
region  and  its  linear  gradient  transition  zone  were  varied  in 
order  to  determine  their  effect  on  the  generation  of  teleseismic 
Rayleigh  waves.  By  comparing  amplitudes  from  a  source  region 
with  sharp  boundaries  at  the  sides  and  bottom  with  those  from  a 
region  bounded  by  various  combinations  of  sharp  and  transition 
boundaries,  it  was  determined  that  for  these  source  dimensions 
at  the  periods  of  interest,  the  more  the  body  wave  energy  that 
leases  the  source  region  as  downgoing  waves,  the  larger  the 
fundamental  mode  Rayleigh  wave. 

In  Section  III,  a  review  on  the  theory  and  application  of 
synthetic  seismograms  is  presented.  The  emphasis  is  on  body 
phase  wave  forms  at  teleseismic,  regional,  and  local  eplcentral 
distances.  At  teleseismic  distances,  it  was  shown  that  long 
period  body  phases  from  shallow  earthquakes  are  coherent  at 
neighboring  stations  and  that  the  observed  waveform  could  be 
decomposed  in  a  manner  that  allows  determination  of  faulting 
parameters.  By  modeling  both  long  and  short  period  body  waveforms 


AFOSR-TR 


o 


^B*5!VrOX£V‘ 


1 


r\ 


using  distributed  point  sources  it  is  found  that  most  earthquakes 
are  Indeed  complex  and  that  high  frequency  strong  motions  appear 
to  be  more  strongly  controlled  by  the  jumps  in  the  source  time 
history  than  by  the  overall  duration.  At  ranges  less  than  30* t 
body  wave  arrivals  become  multi-valued  and  some  care  needs  to  be 
taken  to  remove  propagation  features  from  source  phenomenon. 
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SUMMARY 


The  research  performed  under  the  contract  during  the  period 
1  November  1981  through  30  April  1982  can  be  divided  into  two 


main  topicS/fifcoupling  of  surface  waves  in  laterally  inhomogeneous 
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source  regions  to  teleseismic  propagation  paths,  and  a  review  of 

f 


the  theory  and  application  of  synthetic  seismograms. 

In  Section  II,  the ^dimensions  of  the  cylindrical  source 


region  and  its  linear  gradient  transition  zone  were  varied  in 


O' 


order  to  determine  their  effect  on  generation  of  teleseismic 
Rayleigh  waves.  By  comparing  amplitudes  from  a  source  region 
with  sharp  boundaries  at  the  sides  and  bottom  with  those  from  a 
region  bounded  by  various  combinations  of  sharp  and  transition 
boundaries,  it  was  determined  that, for  these  source  dimensions 
at  the  periods  of  interest,  the  more  the  body  wave  energy  that 
leaves  the  source  region  as  downgoing  waves,  the  larger  the 
fundamental  mode  Rayleigh  wave.  •> 


In  Section  III,  a  review  on  the  theory  and  application  of 
synthetic  seismograms  4s~prese«ted. — The/ emphasis  ^Ls  on  body 


phase  wave  forms  at  teleseismic,  regional,  and  local  epicentral 
distances.  At  teleseismic  distances,  it  was  shown  that  long 
period  body  phases  from  shallow  earthquakes  are  coherent  at 
neighboring  stations  and  that  the  observed  waveform  could  be 
decomposed  in  a  manner  that  allows  determination  of  faulting 
parameters.  By  modeling  both  long  and  short  period  body  waveforms 


using  distributed  point  sources  It  is  found  that  most  earthquakes 
are  indeed  complex  and  that  high  frequency  strong  motions  appear 
to  be  more  strongly  controlled  by  the  jumps  in  the  source  time 
history  than  by  the  overall  duration.  At  ranges  less  than  30  , 
body  wave  arrivals  become  multi-valued  and  some  care  needs  to  be 


taken  to  remove  propagation  features  from  source  phenomenon. 
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II 

Representation  Theorem  Modeling  of  Rayleigh 
Waves  from  Laterally  Inhomogeneous  Source  Regions 

by 

Peter  Glover  and  David  G.  Harkrider 

INTRODUCTION. 

During  the  past  year,  we  have  examined  the  effectiveness  of  two 
mixed-path  surface  wave  algorithms  previously  used  to  couple  various 
source  region  models  to  a  common  propagation  path.  Rayleigh  waves 
generated  by  the  two  algorithms  were  compared  to  the  more  numerically 
accurate  Representation  Theorem  technique  of  modeling  which  couples 
finite  element  code  results  of  a  source  in  a  laterally  Inhomogeneous 
region  to  the  laterally  homogeneous  propagation  path.  The  principal 
advantage  of  the  two  approximations  is  that  they  are  computationally 
economic  and  require  only  trivial  modifications  to  existing  Rayleigh 
wave  codes.  On  thp  other  hand,  the  RT  technique  is  relatively  expensive 
depending  on  the  frequencies  or  wavelengths  desired  from  the  finite 
element/finite  difference  code  modeling  the  transition  cone.  Since 
changing  the  source  region  structure  and  radius  for  each  case  of 
Interest  would  require  further  FE/FD  code  calculations,  we  had  hoped 
that  at  least  one  of  the  two  mixed  path  techniques  would  be  adequate  for 
obtaining  Green's  functions  for  Rayleigh  waves  generated  In  complex 
source  regions. 

In  our  semi-annual  report  for  November,  1981  we  showed  that  for  the 
case  where  the  propagation  medium  was  a  uniform  half space,  the 


Conservation  of  Lateral  Energy  Flux  (CLEF)  method  gave  better  agreement 
with  the  RT  result  than  the  Unit  Transmission  Coefficient  (UTC)  method 
when  the  source  region  was  modeled  by  a  cylinder  of  radius  1.8  km  and 
height  1.8  km  embedded  In  the  propagation  medium.  The  material 
properties  of  this  cylinder  where  chosen  to  closely  approximate  those  of 
two  typical  NTS  source  media.  Climax  Stock  granite  and  Yucca  Flat  tuff 
(see  Table  1).  These  results  are  reproduced  in  Figure  1.  Figure  2 
shows  a  similar  calculation  for  the  case  where  the  propagation  medium 
was  C1T109,  a  model  frequently  used  to  calculate  surface  wave 
propagation  In  the  western  D.S.  Together,  these  figures  show  that  for 
the  Climax  Stock  model,  where  the  source/receiver  rigidity  ratio  is  0.6, 
both  approximations  agree  with  the  RT  results.  However,  for  the  Yucca 
Flat  model,  where  the  rigidity  ratio  is  0.1,  the  CLEF  method  is  the 
better  approximation. 

The  results  given  In  Figures  1  and  2  are  for  a  sharp  boundary.  In 
this  report  we  examine  the  effects  of  gradational  boundaries,  keeping 
essentially  the  same  source  region  geometry  and  using  a  vertical 
point-force  source  at  a  depth  of  0.4  km.  For  the  Yucca  Flat  tuff  source 
region  model,  we  find  that  a  gradational  boundary  on  the  bottom  of  the 
cannlster  gives  the  largest  Rayleigh  waves.  However,  the  maximum 
peak-to-peak  amplitude  Is  6.0  x  10**-21  cm  compared  to  3.5  x  10* *-21  cm 
for  a  similar  force  of  1  dyne  In  the  laterally  homogeneous  model  CIT109. 


EFFECTS  OF  GRADATIONAL  BOUNDARIES. 


Figure  3a  shows  the  RT  and  CLEF  results  for  the  sharp  1.8  x  1.8  km 


Yucca  model  plotted  on  a  common  scale.  The  maximum  peak-to-peak 
amplitude  is  3.9  x  10**-21  cm  for  the  CLEF  method  compared  to  4.0  x 
10**-21  cm  for  the  RT  method.  The  associated  periods  are  15  and  12 
seconds  respectively.  The  CLEF  approximation  underestimates  the 
amplitudes  at  25  second  periods.  At  3  to  5  second  periods,  the  two 
methods  disagree  in  both  amplitude  and  phase.  Nevertheless,  overall  the 
CLEF  approximation  is  in  good  agreement  with  the  RT  result. 

The  underlying  assumption  of  the  CLEF  method  is  that  the  horizontal 
energy  flux  transported  In  the  fundamental  mode  Rayleigh  wave  on  the 
left  of  the  boundary  is  equal  to  that  in  the  fundamental  mode  Rayleigh 
wave  on  the  right  of  the  boundary.  That  is,  there  is  no  energy  loss  due 
to  mode  conversions,  or  reflections.  This  is  not  true  for  the  finite 
element  calculations  of  the  forcing  functions.  IgES  was  dellberately 
chosen  to  be  in  the  propagation  region  in  order  to  monitor  the  change  in 
wave fie Id  between  the  source  region  and  the  the  propagation  medium.  For 
the  high  Impedance  contrast  Yucca  model,  the  outgoing  wavefleld  should 
experience  significant  reflections  and  conversions  on  encountering  the 
sharp  boundary  of  the  cannister.  Therefore,  we  would  expect  the  CLEF 
approximation  to  give  an  even  closer  agreement  with  a  RT  model 
incorporating  an  Impedance  gradient. 

To  investigate  this  effect,  we  Introduced  a  series  of  stepwise 
Increases  in  velocity  and  density  across  the  3  columns/rows  of  elements 
forming  the  RHS/bottom  of  the  plug.  The  radlus/height  was  also 
Increased  to  2.0  km  to  accomodate  this  gradient.  However,  Iggg  was 
maintained  at  r  ■  2.1  km,  z  “  2.1  km  through  these  runs.  Figure  3b 


shows  the  results  for  a  "uniform"  gradient  across  the  RHS  of  the  plug 
compared  to  the  CLEF  result  for  the  sharp  1.8  x  1.8  km  plug.  In  this 
case,  the  RT  results  have  a  peak-to-peak  amplitude  of  3.3  x  10**-21  cm, 
compared  to  3.9  x  10**-21  cm  for  the  CLEF  results.  The  minimum  value  In 
the  CLEF  result  Is  delayed  by  approximately  3  seconds  with  respect  to 
the  RT  result.  Rowever,  the  two  methods  give  nearly  the  same  result  at 
the  longer  periods.  This  is  not  due  to  the  presence  of  the  Jlent 
reducing  the  effective  radius  of  the  sharp  boundary,  as  can  be  st  from 
Figure  3c.  In  fact,  when  we  compare  the  effect  of  plug  radlu  the 
results  using  the  RT  method  alone  (Figure  4a),  It  is  clear  tha  such 
small  perturbations  have  no  measurable  effect.  Therefore,  the  presence 
of  the  gradient  must  be  diminishing  the  signal  amplitudes  by  suppressing 
reflections  from  the  side  of  the  plug.  This  is  further  demonstrated  in 
Figure  4c,  where  the  RT  result  for  the  gradient  model  shows  smaller 
amplitudes  at  all  frequencies  than  the  RT  result  for  the  sharp  boundary, 
and  an  18  percent  decrease  In  the  maximum  peak-to-peak  value. 

Using  a  plug  of  height  2.0  km  and  radius  1.8  km,  we  incorporated 
the  same  velocity  and  density  gradient  across  the  bottom  of  the 
cannlster.  Figure  4b  shows  the  new  result  compared  to  the  RT  method  for 
the  sharp  plug.  The  gradient  gives  larger  amplitudes  at  all  but  the 
shortest  periods  and  the  peak-to-peak  maximum  Is  6.0  x  10**-21  cm  or  50 
percent  larger  than  the  sharp  boundary.  Moreover,  there  is  little 
difference  in  phase  between  the  two  results. 

We  also  computed  the  Rayleigh  wave  for  a  gradient  on  both  the 
bottom  and  the  side  of  the  plug.  These  results  are  shown  in  Figure  5. 
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In  this  case,  the  peak-to-peak  amplitude  is  5.2  x  10**-21  cm,  somewhat 
smaller  than  that  of  the  bottom  only  gradient,  but  still  significantly 
larger  than  for  the  sharp  boundary.  At  other  periods,  the  sharp  vs 
gradational  results  are  in  good  agreement.  Finally,  Figures  5b  and  5c 
compare  the  totally  surrounded  plug  with  the  previous  cases.  The 
Important  thing  to  note  is  that  the  gradient  across  the  bottom  of  the 
plug  has  the  largest  effect  on  the  peak-to-peak  amplitude  of  the 
computed  Rayleigh  wave  signatures. 


DISCUSSION. 


The  observation  that  the  maximum  amplitudes  come  from  the  source 
region  model  with  velocity  gradient  at  depth  can  be  explained  as 
follows.  The  effect  of  the  sharp  boundary  on  the  si  le  of  the  cannister 
is  to  reflect  rays  at  non-normal  angles  of  incidence  so  that  they 
eventually  impinge  on  the  bottom  surface  of  the  cannister.  Here  the 
contributions  due  to  the  normal  stresses  are  a  factor  of  2  larger  than 
those  for  the  side,  this  implies  that,  for  the  periods  of  interest 
here,  the  more  body  wave  energy  that  leaves  the  source  region  as 
downgoing  waves,  the  larger  the  fundamental  mode  Rayleigh  wave.  The 
gradient  itself  acts  so  aB  to  minimize  reflections  from  the  bottom  of 
the  cannister.  The  results  from  the  gradient  on  both  bottom  and  side 
are  consistent  in  that  the  peak  amplitudes  are  somewhat  less  than  the 
case  of  the  gradient  across  the  bottom,  but  greater  than  those  for  the 
sharp  case.  Vhen  the  gradient  occurs  only  on  the  side,  the  results  are 
smaller  than  for  the  sharp  boundary,  which  is  consistent.  Moreover,  the 


V- 


results  for  this  case  most  nearly  match  the  CLEF  results,  although  at 
the  distances  Involved  the  energy  la  being  transported  across  the 
Impedance  contrast  as  body  waves  and  static  deformations,  rather  than  as 
fundamental  mode  surface  waves  (see  Harkrider,  1981  for  a  derivation  of 
the  CLEF  approximation). 

So,  although  our  results  with  these  spatially  limited  source 
regions  clearly  indicate  that,  when  observed  with  a  long-period 
LRSM-type  Instrument,  effects  of  lateral  heterogeneity  can  be  modeled  to 
first  order  using  the  CLEF  approximation,  we  cannot  offer  a  detailed 
explanation  as  to  why. 

REFERENCE 

Harkrider,  D.  G. ,  Coupling  near  source  phenomena  Into  surface  wave 
generation;  Identification  of  Seismic  Sources  -  Earthquake  or 
Underground  Explosion,  ed.  E.  S.  Husebye  and  S.  Mykkeltvelt,  D. 
Reldel  Publishing  Co.,  Boston,  pp.  277-326,  1981. 
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Layer  parameters 

for  the 

Climax  Stock  and 

Yucca  Flat 

models 

a 

6 

P 

P 

Tucca  Flat  Tuff 

2.35 

1.3 

1.86 

3.14 

Climax  Stock  Granite 

5.33 

2.78 

2.67 

20.63 

CIT109  (Upper  14  km) 

6.2 

3.511 

2.736 

33.73 
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CAPTIONS. 

1.  Comparison  of  mixed  path  results  at  1200  km;  source  medium 

1.8  x  1.8  km  right  cylinder,  sharp  boundaries;  propagation 
medium  uniform  half space  a  ■  6.2  km/sec,  8  ■  3.5  km/sec,  p  - 
2.7  gm/cc.  Both  approximations  overestimate  high 
frequencies. 

2.  Similar  calculations  for  layered  propagation  medium.  UTC 

result  for  Yucca  Flat  tuff  still  dominated  by  high 

frequencies. 

3.  From  top  to  bottom,  RT  results  for  a)  sharp  boundaries,  b) 

gradational  boundary  on  side,  c)  sharp  boundaries  but 
reduced  radius,  superimposed  on  CLEF  result  with  maximum  and 
minimum  amplitudes  of  2.0  and  -2.2  x  10**-21  cm. 

4.  RT  results  for  a)  source  region  radius  1.6  km,  b)  gradational 

boundary  across  bottom,  c)  across  side,  superimposed  on  RT 
result  for  sharp  1.8  x  1.8  km  source  region,  peak-to-peak 
amplitude  4.0  x  10**-21  cm. 

5.  RT  results  for  a)  sharp  1.8  x  1.8  km  source  region,  b) 

gradational  boundary  across  side,  c)  across  bottom,  compared 
to  totally  gradational  case,  peak-to-peak  amplitude  5.2  x 
10**-21  cm. 
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THEORY  AND  APPLICATION  OF  SYNTHETIC  SEISMOGRAMS 


Don  V.  Helmberger 


Set  sinological  Laboratory,  252*21 
California  Institute  of  Technology 
Pasadena,  California  91125 


1.  Introduction 

^  ■  ■  .  ■  ■  ■  i  — 

In  the  past  decade,  our  understanding  of  earthquakes  and  earth 
structure  has  Increased  significantly  based  on  our  Improved  ability  to-*- 
Interpret  seismograms.  Numerous  formalisms  have  been  developed  whereby 
synthetic  seismograms  can  be  computed  at  positions  on  the  earth  for 
comparisons  with  observed  seismograms.  The  various  processes  that 
affect  the  motion  at  some  field  point  such  as  the  seismic  source,  earth 
model  response,  attenuation  and  recording  equipment  are  all  combined 
Into  the  synthetic  convolution  operators.  Since  these  effects  can  be 
treated  as  linear  operators  to  first  order  It  becomes  simple  to  test  the 
significance  of  changes  In  the  synthetics  caused  by  varying  the 
operators  separately.  By  comparing  the  synthetic  with  observed 
seismograms  we  can  apply  Iterative  techniques  to  determine  earth 
structure  or  source  model,  or  perhaps  some  properties  of  each  depending 
on  the  circumstances.  Applications  of  this  procedure  to  body  waves  have 
proven  quite  effective  In  determining  the  source  properties  of  shallow 
earthquakes.  He  shall  review  these  studies  In  some  detail  In  this  set 
of  notes. 

There  are  four  basic  operators  which  are  generally  Included  In  the 
generation  of  synthetics.  These  represent  the  seismic  Instrumental 
response,  I(t),  the  attenuation  operator,  A(t),  the  source  operator, 
S(t),  and  a  wave  propagation  operator  or  6reen's  function,  M(t).  The 
synthetic  seismogram,  SS(t),  Is  computed  by  the  triple  convolution 


SS(t)  »  I(t>*  A(t)  «M(t)  *S(t> 


where  the  first  of  these  Is  generally  well  known. 

We  will  be  primarily  Interested  In  events  recorded  by  the  World 
Wide  Seismic  System  Network  (WWSSN).  This  global  network  of  stations  Is 
displayed  In  Figure  1.1.  Each  station  records  three  components  of 
motion,  namely  vertical,  north-south,  and  east-west,  on  photographic 
paper.  Because  the  earth  Is  particularly  noisy  at  about  4  sec,  it  Is 
convenient  to  record  each  component  In  two  pass-bands,  short  period  (SP) 
and  long  period  (IP).  The  gains  or  amplifications  of  (SP's)  are 
normally  between  25  to  200  thousand  while  the  (LP's)  run  at  750  to  3000. 

The  attenuation  operator  does  not  Introduce  any  structure  Into  the 
waveforms  In  most  problems  of  Interest  but  only  smooths  the  results. 

The  operator  which  Is  most  commonly  used  could  be  more  properly  written 
A(t,t*)  where  t*  «  /  ds/Q  where  Q  Is  the  quality  factor  and  the 
Integration  Is  performed  along  the  path  of  the  ray,  see  Carpenter  (1967) 
and  Futtennan  (1962).  Most  recent  estimates  of  t*  are  near  1  for 
P-waves  and  4  or  larger  for  shear  waves  assuming  shallow  events.  The 
effects  of  the  attenuation  operator  convolved  with  the  Instrumental 
responses  are  displayed  In  Figure  1.2.  Note  that  for  a  t*  of  (4)  the 
short  periods  are  greatly  reduced.  Thus,  we  should  not  expect  to  see 
short-period  shear  waves  from  shallow  sources  which  Is  generally  true. 

There  Is  considerable  evidence  for  a  frequency  dependent  attenuation  operator 
particularly  at  the  high  frequences  (say  >  3  hz),  see  Minster  (19781 . 

We  will  be  primarily  Interested  In  the  long-period  seismograms  and  will 
apply  the  A(t*,t)  for  simplicity* 

Calculating  M(t)  can  be  quite  simple  or  complex  depending  on 
circumstances  and  eplcentral  distances.  At  teleselsmlc  distances,  30° 
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Fig  1.2  -  The  first  column  displays 
A(t,t*)  as  a  function  of  t*  assum¬ 
ing  0.5,  1.0,  1.5,  and  4,  from 
top  to  bottom.  The  middle  column 
displays  the  convolution  of 
A(t,t*J  with  the  short  period 
Instrumental  response,  and 
similarly  for  the  long  period 
on  the  right.  The  relative 
amplitudes  are  Indicated  by  the 
numbers  above  each  trace. 


Fig  1.3  -  Example  observations  and  synthetics  at  ranges  where 
the  earth  appears  simple  (>  30  ),  slightly  complicated  (upper- 
mantle  ranges,  30°  to  14°),  and  quite  complicated  (less  than  14°). 


to  90°,  we  can  assume  that  M(t)  Is  a  delta  function  and  the  only 
remaining  complexity  Is  due  to  surface  reflections,  see  Figure  1.3. 
from  about  12°  to  30°,  the  earth’s  upper-mantle  triplications  Introduce 
additional  frequency  dependent  arrivals  making  M(t)  more  Interesting. 

At  smaller  ranges,  the  crust  plays  a  more  Important  role  and  the  guided 
waves  following  the  P-wave  prove  quite  useful.  Methods  for  generating 
M(t)  abound  and  their  descriptions  are  covered  in  great  detail  In  any  of 
the  recent  textbooks.  In  this  set  of  notes,  we  will  apply  generalized 
ray  theory,  6RT,  and  the  Cagnlard-de  Hoop  formalism  since  the  author  Is 
the  most  familiar  with  this  technique  and  It  can  be  used  effectively  at 
all  ranges. 

2.  Source  descriptions  and  generalized  ray  theory 

The  seismic  radiation  field  produced  by  earthquakes  can  be 
represented  by  several  means.  Following  the  stress  relaxation  approach, 
one  assumes  the  initial  stress  and  frictional  conditions  and  performs 
the  proper  dynamics  using  analytical  or  numerical  techniques  to  obtain 
the  displacements,  see  for  example  Madariaga  (1976).  Another 
particularly  useful  approach  Is  due  to  Haskell  (1964),  called  the 
shear-dislocation  model.  Following  this  approach  one  does  not  attempt 
to  understand  the  detailed  mechanics  Involved  In  the  actual  fault  zone 
but  simply  states  that  slip  occurs  on  a  specified  surface,  referred  to 
as  a  dislocation. 

However,  before  getting  Into  the  substantial  complexities  Involved 
In  the  shear  dislocation  formalism  It  appears  useful  to  review: a  much 
simpler  spherically  symmetric  source  or  an  idealized  point  source 


explosion. 

a)  Symmetric  Point  Source 

We  assume  a  homogeneous  fluid  with  cylindrical  coordinates  r  and 
z,  see  Figure  2.1.  This  choice  of  coordinates  will  prove  advantageous 
for  studying  a  layered  earth  In  later  considerations.  We  will,  also, 
take  this  opportunity  to  develop  some  mathematical  tools  which  will 
prove  useful.  The  wave  equation  In  cylindrical  coordinates  with  no 
azimuthal  dependence  is 


Fig  2.1  -  Source-receiver  geometry  (a)  and  complex  (p)  plane  (b) 
with  branch  cut  starting  at  (1/V)  and  running  out  along  the  real 
(p)  *xis. 


(2.1)  d%/dr2  +  1/r  d*/dr  +  d%/d*2  -  l/o2  d2d/dt2  -  0 

where  (+)  Is  the  displacement  potential  and  (a)  the  compressions! 
velocity*  Taking  the  Laplace  transform  of  (2.1)  yields 

*2,2)  d2J/dr2  +  l/r  d?/dr  +  d2?/dx2  -  a2/o2i  -  0 

where  (s)  Is  the  Laplace  transform  variable.  A  solution  of  (2.2)  which 
Is  easily  demonstrated  by  substitution  Is 

?(r,*#»)  •  Je(kr)«±v* 


(2.3) 
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where  (k)  1$  the  horizontal  wave  number,  v  •  (k2  ♦  s2/a2)l/2,  and  J0 
the  zero  order  Bessel  function.  A  more  general  expression  1$ 

m 

(2.4)  ♦<*...•>-*<•>  /  Jo(ki)[*00.'V']« 

o 

where  we  take  the  positive  square  root  of  (v)  and  consider  the  solution 
which  converges  for  large  z.  The  functions  A(k)  and  ?U)  are  arbitrary 
with  respect  to  the  operations  In  (2.2). 

In  the  vicinity  of  the  source  the  displacement  potential  must 
satisfy  the  known  behavior  of  a  spherical  wave  which  is 

0  t  <  R/o 

(2.5)  4a  " 

1/8  f  (t  -  R/«)  t  >  R/o 
Taking  the  Laplace  transform  of  (2.5)  yields 


(2.6) 


♦#(R,a)  -  1/R  r(,/o)R  f  (•) 


The  question  now  is  can  we  find  •  way  to  match  this  condition  (2.6)  with 
the  wave  solution  (2.4)  In  terms  of  4{r,z,s).  This  can  be  accomplished 
by  using  a  trick  originated  by  Lamb  (1904).  The  needed  Identity  Is 


(2.7) 


/"  J.( kr).-v|-h| 


which  Is  easily  proven  with  the  aid  of  the  Fourler-Bessel  theorem.  This 
equation  Is  sometimes  referred  to  as  the  Sommerfeld  Integral 


representation.  Thus,  the  solution  In  cylindrical  coordinates  for  an 
arlbtrary  time  history,  f(s),  becomes 


♦(r,z,«)  -  f(s)  f  Ja(ltr)  idk 

(2.8)  J  °  v 

Of  course,  the  Inverse  Laplace  transform  of  (2.8)  must  yield  the  (2.5) 
description.  We  will  produce  this  result  by  performing  the  so-called 
Cagnlard-de  Hoop  transformation,  a  technique  that  proves  useful  In  the 
layered  earth  model.  We  make  a  change  of  variable  k  *-1sp  which 
transforms  Bessel  functions  to  modified  Bessel  functions,  namely 

(2.9)  J0(-iapr)  -*•  i/x  Ko  (*pr)  -  i/*  Kq  (-*pr) 


The  solution  becomes,  after  Stride  (1959), 


(2.10) 

.  -i  i<i2i 

+lf  6 

f  ro(.pr)  1SE 

where  n  ■ 

(!/•’  -  p!U/2. 

-Iffi 

This  expression  can  be  further  simplified  by  noting  the  synmetry 
across  the  real  (p)  axis  (Schwarz  principle  of  reflection)  obtaining 


If  6 

♦  (r,*,s)  -  f(»)  i(e)  Im  f  K(spr) 

(2.11) 


Equation  2.11  can  be  solved  exactly  but  there  are  some  useful 
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approximations  to  discuss.  Using  the  asymptotic  expansion  (9.7.2 
Abramowltz  and  Stegun)  we  have 


,1/2 


(2.12) 


where  v  •  4n2.  Keeping  the  first  term  only  we  obtain 


(2.13) 


♦  (r.x.s) 


»5T  #-e(pr+n|*-h|)  dj> 
n 


For  simplicity  suppose  we  assume  that  f(s)  ■  1/s  or  that  the 
original  source  time  history  Is  a  step  function.  The  delta  function 
response  can  be  determined  later  by  talcing  a  time  derivative  and  the 
response  for  an  arbitrary  source  by  convolution.  The  type  of  Integral 
occurring  In  (2.13)  occurs  often  in  generalized  ray  theory  and  it 
appears  worthwhile  to  discuss  a  rather  Ingenious  trick  for  obtaining  Its 
solution,  see  de  Hoop  (I960).  From  line  source  theory  we  have 

I* 

(2.14)  +(*,*,•)  -  Im  J  *-*(P*  +  n|*-h|) 

o 

Note  that  this  Integral  looks  much  like  the  formal  definition  of  the 
Laplace  transform,  namely 

S(»)  •  J  a”*1  g(t)dt 

o 

This  Identification  leads  one  to  the  following  change  of  variable* 

(2.15)  t  -  pr  +  (1/a2  -  p2)*  |*-hj 
Solving  for  p(t)  by  applying  the  quadratic  formula  we  obtain 
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p  -fr/R2)t  -  (R2/o2  -  t2)*  |*-h|/R2 

(2.16)  (t  <  R  A») 
n  -  |*-h  j/R2  t  •+  (R2/o2  -  t2)  *  r/R2 

and  . 

p  -(r/R2)t  +  l(t2  -  R2/o2)  %  j*-h|/R2 

(2.17)  (t>R/o) 
n  -  |»-h  |t/R2  -  i(t2  -  R2/o2)  *  r/R2 

Taking  the  derivative  of  (2.16)  we  obtain 


He  now  perform  the  Integration  In  the  *p"  plane  such  that  "t"  will  be 
positive  and  real.  Such  a  path  Is  given  In  Figure  2.1b.  The  contour 
(C)  that  runs,  say  along  the  Imaginary  axis  slightly  to  right  of  the 
origin,  from  (0)  to  (1*  )  Is  deformed  by  analytic  continuation  to  the 
above  contour  r.  The  contribution  from  the  arc  at  (•)  Is  zero  and  there 
are  no  enclosed  poles.  Thus, 


Tm /  +  SE 

c 


(2.20) 
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-  /  •■“"(ilf)  “ 

and  o 

♦  -  I.(±  Jf)  . 

where  H(t)  Indicates  the  step  function. 

Returning  to  (2.13),  with  f(s)  ■  1/s,  we  obtain 


(2.21) 


where 


♦  (r.s.t)  -  /27r“  |  WITT  *  J(t)) 


‘PF-  8) 


and  we  have  used 


S£  (1  //t  )  •  hh% 


The  convolution  operation  Indicated  In  (2.21)  can  be  written 


* 

J  J(t)  dT 


which  can  be  evaluated  either  analytically  or  numerically,  depending  on 
the  complexity  of  J ( t) •  In  this  simple  case,  we  can  use  the 
first-motion  approximation  (Knopoff  and  Gilbert  1959)  and  obtain 

J (t)  -  <i2-R2/a2r*  H(x-R/e) 


•  /i/2  5  (t-R/a)"5* 


Fig  2.2  -  Schematic  diagram  explaining  the  convolution  operation  by 
graphical  means.  We  reverse  the  direction  of  the  square-root  operation 
and  perform  the  multiplication  of  the  two  functions  for  values  of 
t  from  (0)  to  (t)  and  sum.  The  results  are  zero  until  the  two 
functions  overlap. 
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The  solutions  for  more  complicated  models  can  be  easily  obtained 
with  this  technique.  For  Instance,  suppose  we  generate  the  step 
response  appropriate  for  the  Interface  problem  (see  Figure  2.3).  There 
will  be  two  generalized  rays  In  the  upper  half space  required  to  solve 
the  boundaiy  conditions,  see  for  Instance  Spencer  (1960)  or  Helmberger 
(1968).  There  will  be  a  direct  generalized  ray  which  Is  Identical  to 
the  response  just  discussed  plus  a  generalized  ray  that  describes  the 
response  returning  from  the  Interface.  The  latter  response  can  be 
written 

(2.22)  ♦<r,*,t)  -  JTTt  -  U/Jt*  J(t)] 

where 


t  •  pr  ♦  rijCs+B)  ,  n±  •  (“2“  P2) 

J(t)  -  !.(■£-  K(P>  |f) 

2 

R(p)  *  (PjBj  -  +  PjBj) 

and  (dp/dt)  Is  the  derivative  along  the  r  contour,  namely 


nj/(R2/o|  -  t2)1,J  t  <  R/Oj 

i£  . 

dt 

iTjj/U2  -  R2/a2>1/2  t  >  R/aj 


Fig  2.3  -  The  responses  of  J(t)  and  $(t)  as  a  function  of  source- 
receiver  separation.  Note  that  the  graphical  convolution  operator 
at  the  larger  ranges  has  two  parts:  before  t  «  R/o,  the  square-root 
singularities  are  congruent  and  produce  a  log  similarity  at 
t  •  R/ui j  the  second  part  corresponds  to  a  weighted  step  response 
as  in  figure  2.2. 
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where  R  ■  {(z  ♦  h)2  ♦  r2)l/2  which  corresponds  to  the  distance  traveled 
along  a  Snell's  law  reflected  path. 

If  a2  >  there  will  be  a  critical  angle  and  a  head  wave.  In 

this  case  J(t)  will  start  before  t  ■  R/a  ,  since  n2  becomes  complex  at 
p  »  l/a2  and  t£(p  «  l/a2)  becomes 


t 


1  \1/2 

Sj)  (*  +  b). 


Evaluating  (2.22)  numerically,  we  obtain  the  step  responses  displayed  In 
Figure  2.3.  Physically,  this  technique  computes  the  response  by  summing 
over  all  possible  ray  parameters  (p)  that  Interact  with  the  boundary. 

b)  Shear  dislocation  (double-couple)  source 

In  this  section  we  will  examine  some  relatively  simple  dislocation 
models.  Starting  with  Haskell's  representation  for  shear  faulting  It  Is 
relatively  easy  to  devise  the  displacements  for  double  couples  In  an 
Infinite  medium,  see  Harkrlder  (1976).  The  solution  for  a  strike-slip 
fault  becomes 


Fig  2.4  -  Coordinate  system  for  the 
dislocation  formulation  where  (e) 
Indicates  the  strike,  («)  the  dip 
of  the  hanging  wall,  and  (x)  the 
motion  of  the  hanging  wall  relative 
to  the  footwall. 


(2.23) 


£)]  .in  29 

v(r,*,«)  -  k  [f-  (0  -  7  If)  +  k|  If]  cos  26 
Q(r.*,«)  -  K  (0  -  F  If)  +  k9  If]  tin  26 


where  N,  V  and  Q  are  the  displacements  In  the  vertical,  azimuthal  and 
radial  coordinates.  The  parameters  are 


▲ 


•  it 

•  s 


B 


& 


where  R  ■  distance  between  the  source  and  receiver;  a, 6  *  congressional 
and  shear  velocities,  and  K  ■  source  strength.  The  spatial  derivatives 
are  Introduced  because  of  the  double-couple  nature  of  the  source,  and 
the  A  and  B  are  simple  point  source  representations  as  before.  Using 
the  Sommerfeld  Integral,  (2.7),  and  performing  the  appropriate  algebra 
(see  Helnberger  1974)  we  can  express  the  displacements  along  the 
vertical,  tangential,  and  radial  coordinates  as: 


(2.24) 


V  -  ||>  +  »pO 


I  11 .  _i_  _  li 

F  39  apr  9x99  9r 

9£_  i_  a2n  l 

9r  "  «p  9r9x  r  99 


Q 


o» 


where  2,  r  and  0  are  the  vertical,  radial,  polar  angle  coordinates 
respectively.  The  P  wave  potential  (♦'),  the  SY  wave  potential  (B),  and 
the  SH  wave  potential  (x)  are  expressed  by; 
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P-wave; 


Mo 

Z7F 

2 

7 

I« 

-H-+c 

/ 

Cjtp) 

c 

i 

4*p 

2 

"* 

Tm 

■H-+C 

/ 

C2(P) 

c 

Mo 

4*9 

2 

* 

Im 

+i»+c 

/ 

c3(p) 

SY-waves: 


-H-+c 


+  W  7  I"  /  w,(P)  ^ 


+1-+C 


7  I"  /  tv,cp)  i 


*I7  7h  /  *»»»  ^ 


SH-waves; 


+1-4C 


i  + zt!  f  *■  / 


+i»fc 


Jm 

4*p  * 


/  “a**1  ^ 


r-hljK^spiOdp'A^e.X.fi) 

*-h  ( )  Kj  (spr )  dp  •  Aj  (6 ,  X ,  6) 

*-h|)Ko(8pr)dp*A3(6,X,6) 

Jr-hll^CsprJdp-Ajfe.X.d) 

r-hllKjCsprldp'Aj 

|x-h))K0(Bpr)dp*A3 

-«nBl*-hjK2*pr)dp*A4 

-sngU-hllK^BprJdp-Aj 


(2.25) 


s  ■  Laplace  transform  variable 
p  ■  ray  parameter 
\  .  (i/v2  -  p2)l/2 
h  «  depth  of  source 
o  ■  congressional  velocity 
8  ■  shear  velocity 
p  •  density 

%  *  seismic  moment 

with  the  orientation  constants  given  by: 

*l(e,x,6)  ■  sin  2e  cos  x  sin  6  «•  1/2  cos  2e  sin  x  sin  26 
^(e.x.d)  ■  cos  e  cos  x  cos  6  -  sin  e  sin  x  cos  2  6 
A3 (0  »x ,  6)  ■  1/2  sin  x  sin  26 

*4(8, x,i)  «  cos  28  cos  X  sin  6  -  1/2  sin  28  sin  x  sin  26 
*5(8, x, 6)  ■  -sin  e  cos  x  cos  6  -  cos  8  sin  x  cos  26 

(2.26) 

where 

•  •  strike  from  the  end  of  the  fault  plane 
X  »  rake  angle 
6  ■  dip  angle 

The  vertical  radiation  patterns,  as  will  become  apparent  shortly 
are  defined  by 


c,  - 


sv,  •  -epn. 


sh,  - 


i  F 


(2.27) 


C2  “  2cpna 


c  ne 


SV2  -  (n l  -  p!>  SH2  -  Jj  f 


C3  -  (p2  -  2r>2)  SV3  -  3cpne 


where  t  ■  { 


+1  *>h 


-1  *<h 


The  displacement  field  for  any  orientation  of  a  double  couple,  see 
Figure  2.4,  Is  the  sum  of  displacement  fields  of  three  particular 
orientations  of  the  double  couple,  corresponding  to  the  strike-slip 
fault  (X  -  0  or  180°,  «  .  go°)#  normal  dip-slip  fault  (90°,  90°),  and 
the  dip-slip  fault  with  a  45°  dip  (45°,  90°),  see  Burrldge  filaL, 
(1964). 

The  Index  (n)  In  the  expression  (2.27)  corresponds  to 
(1)  strike-slip,  (2)  dip-slip,  and  (3)  45°  dip-slip. 

The  Integrals  expressed  In  (2.25)  can  be  transformed  back  Into  the 
time  domain  by  the  application  of  the  Cagnlard-de  Hoop  technique,  see 
Gilbert  and  Helmberger  (1972)  and  Harkrlder  and  Helmberger  (1977).  For 
example,  the  field  function  defined  by 


(2.28) 


tn(r,z,») 


•  <|  a  Im  J  -jE-  Kn(spr)  e“8nvl*"hl  dp 


becomes 


(2.29) 


2  9  -  f  Cn(t,T)  /d£\ 

J  w 


Cn(r,x,t) 


\ 


Vv-' 

CilUk  ii  -  • 


' •'  i*.  **-  1*.  <*  •. 


where 


cn{t,T(p))  -  cosh  cosh"1  ^■t~T^pr 

The  geometry  Is  given  In  Figure  2.1a  and  the  de  Hoop  contour  In 
Figure  2.1b,  see  de  Hoop  (1960).  The  various  functions  of  p  are  to  be 
evaluated  along  r  defined  by  choosing  those  values  of  p  which  make  (t) 
real  and  Increasing  where 

T(p)  -  pr  +  njs-fcl* 

The  transformation  of  "p"  to  "t"  follows  from  the  algebra 
discussed  In  the  last  section. 

In  this  simple  case  we  have  a  closed  form  solution  for  various 
values  of  n  since  the  equivalent  form  back  In  the  («,k)  domain  has  been 
evaluted  by  Harkrlder  (1976).  For  example,  v 

tjCr...*)  -  £  j  [I  +  f|  (t -!)]*(« -|)| 

where  the  near-field  contribution  appears  In  terms  of  r.  However,  since 
we  need  to  evaluate  Integrals  similar  to  (2.25)  with  complicated  complex 
Integrands  later  It  should  be  noted  that  (2.29)  can  be  evaluated  for 
various  values  of  (t)  after  a  change  of  variable  as  discussed  earlier. 

A  relatively  fast  evaluation  of  this  type  of  Integral  Is  by  nonuniform 
quadrature  techniques  where  the  point  spacing  Is  determined  by  the  rate 
of  change  of  the  Integrand.  The  accuracy  of  such  techniques  are 
discussed  In  Helmberger  and  Harkrlder  (1978).  He  will  examine  some 


useful  approximations  by  expanding  the  integrand  of  (2.29)  in  terms  of 
(t  -  t)'\  Note  that 

£  ,t  T  -  1  ffy+(y2-D1/2l2n  +  l 

n  2  L  (y  +  (y2-D1/2)  J 

where 

y  -  (t  -  t  +  pr)/pr 

and  # 

cn(t,T.p)  ! 

1/2  m 

<t-T  +  2pr)A/A  *5pr 


to  first  order.  Thus,  we  can  approximate  (2.29)  by 

t 


-Ul.  f  — zr - 4 


_  J2L  dT 
^2  dt  nvdT 


(2.30) 


£  [f  *  ^  *  %  &)] 


This  expression  was  obtained  earlier  by  using  the  first  term  of 
the  asymptotic  form  of  the  modified  Bessel  function.  With  this  degree 
of  accuracy,  called  the  high-frequency  solution,  the  expressions  (2.25) 
can  be  greatly  simplified.  We  obtain  v 


4  "  4“  Aj(8,X,6)  ^  Im  J  cj  ^  A*/2spr)  *  lexp(-8(pr+njz-h|)dp] 

J-l  o  a 


M  3 

j-1  O  B 


\iZ2spt)  x  [exp(-s(pr+nA|z-h|))dp] 


M  ^ 

X  "  Tip  A^+3)(e,X,6)  |  Im  J  SHj  -2-  /(it/2spr)  *[exp(-s(pr+np|z-h|)) 
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where  we  are  essentially  assuming  that  the  source  duration,  7,  Is 


T  <  <  2  p  r. 


We  still  further  approximate  (2.30)  by  assuming 


dp/dt  - 


<t-tR)*(2tR) 


where  t^  *  R/V  and  p  *  pQ  and  (2.30)  reduces  to 


Cn(r,*,t)  -  6(t-t_)/R 


called  the  first  motion  approximation.  This  approximation  Is  valid  at 
teleselsmlc  distances  where  the  ratio  of  travel  time  to  duration  Is  of 
the  order  of  100  or  greater  and  has  proven  quite  useful  In  modeling 
shallow  earthquakes,  see  Langston  and  Helmberger  (1975). 

Transforming  the  expressions  (2.31)  Into  the  time  domain  we 


obtain: 


♦  -w  Lve-M) 


(2.32) 


3-1 


*  -  Z  y3«.M>s»i 

3-1 


ft  -  (R/< 

R 


Thus,  the  potentials  contain  the  classical  vertical  radiation  patterns 
in  terns  of  the  C^(p),  SVj(p),  and  SH^(p)  and  the  horizontal  pattern  In 
terms  of  the  Aj's.  We  will  use  these  expressions  In  discussing  the 
teleselsmlc  results  later. 

c)  Response  of  a  layer/half space:  Digression 
Using  the  concepts  of  generalized  ray  theory  and  retaining  only 
the  first  term  of  the  asymptotic  solution,  we  can  construct  the  (SH) 
solution  on  a  free  surface,  namely 

(2.33)  V(r.O.«,t>  £  i(t>  *  £  Aj+3<0.M>  Vj(t>  . 

0  L  j-1 

where 

yt)  -  2  Im  (  ^  SHjIplnCpip  J 

D(t)  -  dislocation  history 
D(t)  »  far-fleld  time  function- 

nit)  ■  product  of  the  transmission  and  reflection  coefficients 
Z  ■  summation  over  contributing  rays 

For  the  case  of  a  pure  strike-slip  dislocation  embedded  In  a  homogeneous 
half space 


(2.34) 


Vj  -  Vx  -  (2  •ini/r)H(t  -  tR)/R 


where  (1)  Is  the  angle  of  Incidence,  see  Figure  (2.1a)  .  Substituting  (2.34) 
Into  (2.33)  we  obtain 


(2.35)  V  -  (Mo/4irp  )(2F  8±ni/B^)cos  26  D(t  -  tB)/YR 

-20 

where  the  amplitude  Is  given  In  centimeters  with  F0*  10  for  unit 
conversion,  and  the  various  parameters  expressed  as  M0(ergs)(  e0(g/Cjn3)„ 
e(lcm/sec),  R(km)  and 

Y  mJ‘  D(t)dt  -  »(•) 

A  factor  of  2  has  been  Introduced  In  the  solution  (2.35)  by  the  SH 
free  surface  receiver  coefficient.  The  far-fleld  step  function  response 
for  a  pure  strike-slip  event,  Yj(t),  is  given  In  Figure  2.5,  where  the 
model  Is  Included  as  an  Inset  for  various  values  of  source  depth.  The 
corresponding  half-space  response  Is  the  simple  step  displayed  on  each 
trace.  Comparing  the  response  with  the  source  situated  just  above  (d  « 
3.5)  and  below  (d  «  4.5),  one  finds  that  the  long-period  behavior  Is 
nearly  the  same.  The  classical  type  of  Love  wave  dispersion  Is 
developed  when  the  layer  contains  the  source  and  Is  well  understood  In 
terms  of  ray  Interference.  When  the  source  Is  located  below  the  layer, 
the  Interpretation  Is  more  difficult  but  can  be  studied  by  examining  the 
various  rays.  The  ray  descriptions  become: 


pr  +  n2 

(d-h)  +  n2h 

pr  +  n2 

(d-h)  +  nx(3h) 

pr  +  n2 

(d-h)  +  nx(2n-l)h 

T21<»> 

In<P)  ! 

«12(P> 

T2i<»  ' 

d«4.5  krrt 


0*2J5  p* 2j6  h»4  km 


U 

Source  depth  «d,  A*7Z  km 


Time,  sec 

Fig  2.5  -  SH  step  responses  at  the  surface 
assuming  a  point  strike-slip  dislocation 
situated  at  various  depths.  The  amplitudes 
are  scaled  relative  to  the  top  trace  with  the 
step  response  for  a  homogeneous  half-space 
(bottom  properties)  displayed  for  comparison. 
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0*  (p) 

Fig  2.6  -  Contours  of  I m  versus  Re(p)  for  the  direct  rays 
shown  in  Figure  2.5.  The  branch  cuts  run  along  the  Re(p) 
coordinate,  starting  at  l/e2  and  1/Bj,  respectively.  The 
parameter  time  Is  marked  along  each  contour. 

where  the  index  (1)  refers  to  the  layer  parameters  and  (2)  to  the  bottom 
half  space. 

The  coefficients  are: 


end 

where 


1^2  (p)  -  (v^  -  w2n2)/(y1n1+  v2n2) 
T21(p)  -  (2w2n2>/(Vi  +  »»2n2) 


wl  "  P1®1  • 


w2  -  P2e2 


Each  ray  In  the  solution  must  be  evaluated  along  Its  own  contour 
which  Is  determined  by  Inverting  these  polynomial  equations, 

P(*n)»where  t  is  positive  real  and  Increasing  from  0  to  ••  The  contours 
used  to  compute  the  direct  rqys,  (ty,  for  the  responses  In  Figure  2.5 
are  shown  In  Figure  2.6.  The  parameter  (t)  Is,  also,  plotted  along 
these  contours  and  It  Is  easy  to  see  that  abrupt  Increases  In  Im(dp/dt) 
lead  to  arrivals  In  Figure  2.5. 

Note  that  lm(dp/dt)  always  starts  with  the  square-root 

singularity.  This  singularity  yields  the  geometric  ray  arrival.  For 

large  source  depths,  the  contour  Is  near  vertical  and  the  synthetic 

waveform  closely  resembles  the  geometric  rty  response.  The  first  motion 

approximation  Is  actually  equivalent  to  the  saddle  point  approximation 

which  assumes  that  the  contour  goes  straight  up  to  (-).  Thus,  the  first 

motion  approximation  Is  useful  for  teleselsmlc  studies.  An  example 

application  of  this  technique  of  earthquake  modeling  Is  given  In 

Figure  2.7.  The  observed  displacement  Is  from  a  simple  strike-slip 

earthquake  occurring  In  the  Imperial  Valley,  California.  The  only 

unknown  parameters  in  this  particular  situation  were  the  depth  and  slip 

history,  D(t).  After  a  diligent  search,  the  source  depth  of  7  km  and  a 
• 

D(t)  specified  by  triangular  pulse  with  duration  1.5  sec  was  determined 
(see  Heaton  and  Helmberger  1978).  Similar  studies  have  been  conducted 
by  Heaton  and  Helmberger  (1977)  and  Helmberger  and  Malone  (1975).  A 
rather  common  observation  In  modeling  multi-bounce  phases  Is  that  they 
tend  to  change  their  character  relatively  rapidly.  For  Instance,  many 
times  the  first  few  bounces  will  appear  nearly  Identical  at  neighboring 


Fig  2.7  -  Comparison  of  a  synthetic  with  a  strong 
motion  recording,  IVC,  of  an  earthquake  occurring 
at  Brawley,  California,  In  November  1976.  The 
model  parameters  expressing  the  known  structure 
are  given  at  the  top  of  the  figure  with  the  source 
at  a  depth  of  6.9  km.  The  three  most  Important 
rays  are  Indicated  schematically. 


stations  as  predicted  for  flat  layers.  However,  the  next  bounce  will 
suddenly  be  missing  or  too  large  at  these  same  stations.  Another  common 
feature  Is  a  sudden  change  In  frequency  content  of  the  multiples.  Such 
effects  could  be  caused  by  lateral  variation  In  the  low- velocity 
waveguide  near  the  surface.  A  simulation  of  the  effects  of  a  simple 
dipping  Interface  Is  displayed  In  Figure  2.8  (see  Hong  and  Helmberger 
1977),  In  this  particular  case,  the  rays  can  go  beyond  critical  angle 
after  bouncing  and  thus,  produce  high  frequency  reflections.  Some 
progress  on  treating  smoothly  varying  3-dlmenslonal  structure  is 
presented  by  Hong  and  Helmberger  (1977). 

d)  Full  Cagnlard  Solution 

The  high  frequency  solutions  given  by  (2.31)  have  many  advantages 
In  model  studies  as  just  discussed.  However,  for  small  values  of  (spr) 
one  must  use  the  full  solution  by  applying  the  transformation  given  in 
expression  (2.29).  The  displacements  are  obtained  by  substituting 
expressions  (2.25)  Into  (2.24)  and  performing  the  Cagnlard-de  Hoop 
technique.  The  results  are  complicated  because  of  the  near-field  terms 
The  tangential  motions  become: 

v(i,i.e,t)  -  li(t)*  pf  *3+3  vj ’ 

Vi(r,«,S,t)  -  4  X"  J  «e<2)  SBjBj  it 

♦  \imir  «s(2)  dT  4* 

♦  4  'mfT  «.<»  clRFT(r)  *  * 

♦  \uff  *.«>  sVst<4>  4t  dt 
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Fig  2.8  -  Comparison  of  the  step  responses  for  a 
flat  and  dipping  layer  where  the  elastic  parameters 
are  the  same. 


Asymptotic  approximation 


47 


Point  strike -slip  dislocation 


Tangential 


Vertical 

^ 


t-S 


Fig  2.9  -  Comparison  of  the  various  components  of 
motion  for  a  strike-slip  orientation  at  A  «  16  km 
for  a  vhole  space.  The  top  four  rows  contain  the 
asymptotic  sunwatlon  after  1,2,4,  and  12  Items.  The, 
full  solution  Is  displayed  on  the  bottom. 


where 

-  p,  Rpj  -  1,  *ST  -  -ng/p  for  whole  epece 
and 

4n  n»  "2na  <*!  -  ?2> 

*--2p.  *?t  --r-*’  rst  “  "T^  — 2 - 

T  **  62R(p)  ST  p  6  R(p) 


for  a  receiver  on  the  free  surface.  The  dip-slip  result  Is  similar  with 
n  *  1  and  where  the  factor  In  parentheses  Is  reduced  by  two. 


An  evaluation  of  this  expression  for  a  whole  space  Is  given  In 
figure  2.9,  where  the  near-field  contributions  are  apparent.  The 
asymptotic  solution  Is  Included  for  comparison.  These  computations  are 
expensive  and  their  applications  to  data  complicated,  see  for  Instance 
Helmberger  and  Harkrlder  (1978)  and  Heaton  and  Helmberger  (1979).  Thus, 
we  will  limit  our  discussion  to  the  more  useful  approximations  where 
near-field  terms  can  be  neglected. 

3.  Modeling  at  Teleselsmlc  Distances 

Synthetic  seismograms  can  be  constructed  for  comparison  with 
teleselsmlc  observations  by  putting  together  the  various  sub-operations 
discussed  In  the  previous  sections.  Ne  will  be  primarily  Interested  In 
shallow  events  where  Independent  Information  Is  many  times  available. 
Thus,  we  must  correct  for  the  free  surface.  The  simplest  correction  for 
this  situation  Is  to  assume  Ideal  elastic  Interaction  and  Include  the 
reflected  phases,  namely  pP,  and  sP  (see  Figure  3.1).  Using  the 
expressions  (2.24)  and  (2.32)  and  the  concepts  of  GRT  (see  Appendix  A) 


Fig  3.1  -  Schematic  representation 
of  the  phases  P,pP,  and  sP. 


/ 

i 


we  can  approximate  the  vertical  response  at  the  receiver. 


(3.1) 


0  j-i  L 

n_  1  r|pz|  <P0> 

+  *W  T&  R.p(P0)6(t-At25J  Aj  — B— 


where 

atA  ■  2Hna,  time  lag  of  pP 
tt2  ■  H(na  +  tig),  time  lag  of  sP 
H  -  source  depth 
D(t)  ■  6(t).  delta  function 

The  expression  In  brackets  can  be  represented  by  three  spikes  with 
various  strengths,  depending  on  the  source  orientation,  referred  to  as 
stick  diagrams,  see  Figure  3.2.  The  diagram  given  In  the  first  column 
displays  this  Information,  along  with  the  timing  of  the  three 
Interacting  phases,  P,  pP,  and  sP,  for  the  three  fundamental  faults. 
Crustal  layering  Is  easily  Incorporated  by  applying  ray  summations  as 
discussed  In  the  previous  section.  We  assume  •  1  and  that  the  (1/R) 
geometric  spreading  term  In  the  homogeneous  case  can  be  replaced  by  the 
effective  •1/R"  from  Figure  3.3.  The  synthetics  In  Figure  3.2  were 
produced  by  performing  the  various  convolutions  assuming  different 
source  durations. 

The  waveshapes  are  seriously  distorted  by  the  interferences  with 
each  fault  orientation  having  Its  own  characteristics.  Thus,  we  can 
compare  these  synthetics  with  observed  waveshapes  such  as  those  given  In 
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(M0/4wp  )-  (1/R),  and  receiver  functions  Included.  The  source 
time  parameters,  St's,  are  high  stress-drop  (0.5,  1.0,  0.5), 
medium  stress-drop  (1.0,  3.0,  1.0),  and  low  stress-drop 
(2.0,  6.0,  2.0).  After  Langston  and  Helmberger  (1975). 


51 


Figures  3.4,  3.5  and  3.6  to  determine  the  fault  parameters;  namely, 
strike,  dip,  slip  vector,  fault  depth,  moment,  and  fit's  (source  time 
history).  Note  that  for  observations  near  nodes  such  as  BOG  and  BHP 
(Figure  3.4)  one  finds  difficulty  In  determining  the  polarity  of  the 
direct  P  since  It  Is  nearly  zero.  In  cases  where  fi  f  90°  we  expect  sP 
to  still  contribute  as  at  BOG  since  the  dip-slip  component  Is  strong. 
Note  the  strong  differences  In  amplitudes  displayed  In  Figure  3.2  which 
are  caused  by  the  position  of  the  radiation  loops. 

Considerable  ringing  occurs  for  those  events  situated  In  some 
basins  and  In  oceanic  environments  (see  Figure  3.6  and  Langston  and 
Blum,  1977,  for  examples  of  multl-pathlng  or  situations  where  many  rays 
are  Involved  in  generating  the  response). 

Synthetics  for  normal  and  thrust  faults  are  not  very  Interesting 
since  they  are  nearly  Identical  at  all  azimuths.  As  probably  suspected. 


Borrego  Ml 


Fig  3.4  -  Observed  (top)  and 
synthetic  (bottom)  long  period  P 
waveforms  at  14  WMSS  stations. 

The  P-first  motion  plot  Is  rep¬ 
resented  by  the  equal  area 
stereographic  projection  of  the 
lower  half  of  the  focal  sphere. 

Black  dots  Indicate  compression 
(upward  breaking  P)  and  open  circles 
indicate  dilatation  (downward 
breaking  P).  The  heavy  solid  lines 
denote  the  nodal  planes  used  In 
determining  the  fault  orientation, 
e  (strike),  6(dip),  x(slip  direction). 
Modified  after  Burdick  and  Mellman 
(1976). 
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Fig  3.5  -  Focal  mechanism 
plot  showing  the  P-flrst 
motions  for  the  Orovllle 
earthquake,  August  1,  1975. 
The  observed  P  waveforms 
are  given  on  the  top  and 
synthetic  results  on  the 
bottom.  Modified  after 
Langston  and  Butler  (1976). 
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this  uniformity  In  waveshape  leads  to  a  lack  of  resolution  In  fault 
orientation,  ire  require  samples  of  the  energy  release  at  other 
quadrants.  In  particular  at  regional  distances.  At  these  distances  the 
earth's  mantle  structure  produces  more  complicated  responses  and  the 

1/R  6(t~tg)  correction  Is  no  longer  appropriate. 


4.  Modeling  Regional  Bo4y-waves 

The  regional  phases  can  be  subdivided  Into  two  groups,  those  less 
than  and  those  beyond  about  12°,  At  the  smaller  ranges  the  response  Is 
controlled  largely  by  the  crustal  waveguide  and  at  larger  ranges  by  the 
upper-mantle  discontinuities,  see  Helmberger  and  Engen  (1980). 


54 


a)  Crustal  waveguide 

We  begin  by  discussing  the  effects  of  the  crustal  waveguide  on  the 

P-waves,  actually  all  the  energy  that  travels  faster  than  the  S-waves.  This 

portion  of  the  seismogram  Is  conmonly  referred  to  as  PR,  PL,  and  P  ..  We 

n  pul 

assume  a  point  shear  dislocation,  as  before,  and  that  the  earth  can  be 
treated  by  a  single  layer  (crust)  over  a  half-space  (mantle),  see  Table  1. 
Both  of  these  assumptions  are  obvious  oversimplifications  but  as  we 
will  see,  worth  consideration  as  viable  models  at  long  periods.  The 
vertical  displacement,  assuming  the  high-frequency  solution,  becomes 


(4.1) 

where 


(4.2) 


Vj(t) 


tP|  *  L  tm  (  \  C1«  &)] 


where  R^fp)  Indicates  the  appropriate  receiver  function  for  either 
P(N*P)  or  S(N«S)  waves  arriving  at  the  station,  respectively.  The  radial 
displacements,  Qj,  are  obtained  by  replacing  and  RpR  with  R^  and  RSR, 
see  Appendix  A  for  definitions  of  these  receiver  functions.  The 
sumnatlon  Is  over  contributing  rays  where  n^(p)  defines  the  product 
of  all  the  reflection  coefficients  encountered  along  the  ray  path  from 
the  source  to  the  receiver,  similar  to  the  SH-case  considered  earlier. 


Constructing  the  solution  by  sunning  the  rays  as  described  by  (4.2)  is 
tedious  in  that  many  rqys  are  required  to  produce  reliable  synthetics. 

To  ensure  stability*  ray  files  containing  increasing  numbers  of 
multiples  are  computed  sequentially  where  the  output  can  be  monitored 
per  multiple.  An  example  calculation  is  displayed  in  Figure  4.1.  Since 
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Fig  4.1  -  Vertical  components  of  motion  as  a  function 
of  ray  summation  assuming  a  trapezoidal  time  function 
described  by  0.5,  0.5,  0.5,  for  pure  strike-slip  and 
dip-slip  orientations  at  a  ■  1000  Ion. 
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the  strengths  and  polarities  of  the  various  rays  are  strongly  Influenced 
by  the  vertical  radiation  patterns,  Cj(p)  and  SVj(p),  we  have  Included 
these  patterns  above  each  column.  Note  that  due  to  our  limited  time 
window  we  do  not  Include  the  direct  S  pulse.  The  second  pulse  arriving 
In  the  direct  ray  response  Is  actually  a  diffracted  arrival,  essentially 
the  pulse  that  travels  to  the  surfaces  as  an  S-wave  and  along  the 
surface  as  a  P-head  wave.  The  first  set  of  multi-bounce  rays  produce 
head  waves  along  the  top  of  the  mantle.  There  are  12  rays  In  this  set, 
namely 

22  ♦  23  .  12 

considering  the  various  possible  mode  changes.  The  responses  In  the 
bottom  row  were  generated  by  summing  all  (4,094)  rays  as  a  check  against 
a  much  smaller  ray  set  which  takes  advantage  of  the  kinematic  and  dynamic 
redundancies.  For  Instance,  only  304  rays  are  needed  to  compute  the 
response  for  the  fifth  bounce  since  most  of  the  rays  In  the  complete  set 
have  the  same  times  and  amplitudes  (see  Hron  and  Kanasewlch  1971  or 
Helnterger  1968  for  details). 

Another  strategy  for  understanding  these  motions  Is  to  split  these 
ray  contributions  Into  those  starting  upward  as  opposed  to  downward  (see 

Figure  4.2).  Changing  the  source  depth  shifts  the  upgolng  trace 
relative  to  the  downgoing  trace  producing  a  dramatic  effect  on  the 
dip-slip  case  as  displayed  In  Figure  4.3.  In  fact,  a  good  approximate 
seismogram  for  h  •  4  and  16  can  be  obtained  by  simply  shifting  the  h  •  8 
traces  without  recomputing  the  rays.  Partial  derivatives  with  respect 


STRIKE-SLIP  DIP- SLIP 


H30mcH  '  ^ 


Fig  4.2  -  Radial  components  of  motion  at  L  *  1000  tan  in  terms 
of  upgolng,  downgoing,  and  sunned  for  trapezoidal  time 
history  (0.5,  0.5,  0.5  sec). 
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Fig  4.3  -  Motions  of  all  three  types  of  faults  at  A  »  1000  tan 
as  a  function  of  source  depth.  The  numbers  above  each  trace 
Indicate  the  zero  to  peak  amplitude  In  centimeters  assuming 
M0«  4itp.  x  1020  dyne-cm  and  trapezoidal  time  history 
(0.5,  0?5,  0.5  sec). 


to  changes  In  velocities  or  changes  In  the  waveguide  can  be  approximated 
by  stretching  or  compressing  Individual  ray  responses  (see  Wallace  and 
Helmberger  1982). 

For  the  purpose  of  studying  source  mechanisms  we  can  apparently  do 
reasonably  well  by  assuming  that  common  shallow  earthquakes  occur  at  a 
depth  of  8  km.  Profiles  of  responses  for  a  dislocation  located  at  this 
depth  are  presented  In  Figure  4.4.  Preliminary  comparisons  of 
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Fig  4.4  -  Profiles  of  the  vertical  displacements  for  the  three 
fundamental  faults.  The  numbers  above  each  trace  Indicate  the 
zero  to  peak  amplitude  In  centimeters  assuming  Mg  ■  4irp0  x  1020 
dyne-cm  and  trapezoidal  time  history  (0.5,  0.5,  0.5  sec).  This 
time  history  was  convolved  with  the  delta  function  response  to 
reduce  the  high-frequency  spikes  created  by  the  post-critical 
reflections  so  that  the  longer  periods  become  more  apparent 
(After  Wallace  et  al_. ,  1981a). 


synthetics  generated  from  these  responses  (see  Wallace  et  al.,  1981), 
Indicate  that  these  results  are  useful  on  a  global  scale  If  the  high 
frequencies  are  removed.  Thus,  we  will  filter  these  responses  and  the 
corresponding  observations  by  convolving  with  a  triangle  which  has  a 
2  sec  rise  and  fall.  Convolving  these  filtered  Green's  functions  with  a 
trapezoid  source  with  1  sec  Intervals  followed  by  the  long-period  WWSSN 
Instrumental  response  yields  the  synthetics  given  In  Figure  4.5.  This 
time  history  Is  expected  for  events  In  the  5  to  6  magnitude  class,  the 
size  commonly  recorded  by  the  network.  We  propose  to  model  shallow 
earthquakes  at  these  ranges  by  simply  forming  a  linear  combination  of 
these  three  fundamental  synthetics.  The  relative  strengths  of  each 
column  Is  used  to  fix  the  fault  orientation,  namely  the  strike,  dip  and 
slip  vector.  Some  example  comparisons  of  observations  and  corresponding 
synthetics  are  displayed  In  Figures  4.6,  4.7  and  4.8.  These  results 


Fig  4.5  -  Theoretical  displace¬ 
ment  profiles  for  the  vertical 
component.  The  Green's  functions 
were  computed  from  the  model 
presented  in  Table  1  and  have 
been  convolved  with  a  source 
time  function  represented  bv  a 
trapezoid  (ti»l,  t2*l.  t3*l). 
a  triangular  filter  (2  second 
rise  and  fall),  and  a  WWSSN 
long-period  Instrument. 


ZSS  ZDS  Z45 

vfy'' 

Aft  5.4 5.4  A/V 

VYV  ^  63-tyV 

VY\/  72  72  -\V\ 

V'yv  m  ai  Vv^ 

-/yVv  a9-^|V^  as-fy A\f^ 

87  yw 

■/\fAy\l2e  IZe-VVA/ 


Truckee  earthquake 


Fig  4.6  -  The  vertical  waveforms  of  the  1966  Truckee  earthquake. 
The  star  denotes  the  epicenter.  The  data  Is  the  top  trace  at  each 
station  and  the  trace  below  is  the  synthetic  fit.  The  strike-slip 
mechanism  has  two  nodal  planes  which  project  through  TUC  and  BOZ. 

To  the  right  of  each  trace  Is  the  observed  or  predicted  amplitude 
(on  the  basis  of  M  ■  .8  x  1025  dyne-cm)  in  10"3  an. 


Fig  4.7  -  Filtered  data  and 
synthetics  from  the  Orovllle 
earthquake.  At  all  the 
stations  except  GOL  both 
the  vertical  (the  first 
trace  pair)  and  radial 
components  are  shown. 


Tig  4,8  -  Location  of  the 
Turkey  event  (star)  and  the 
recording  stations.  The 
filtered  data  and  synthetics 
for  both  the  vertical  and 
radial  components  are  shown 
(After  Wallace  and  Helmberger, 
1982). 
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were  found  by  applying  a  least  squares  waveform  Inversion  technique. 

The  method  Involves  an  Iterative  scheme  which  makes  use  of  an  error 
function  determined  by  the  cross-correlation  of  the  observed  seismogram 
and  a  synthetic.  Since  only  three  parameters  are  Involved,  the  scheme 
proves  very  effective.  The  moments  are  determined  by  simply  comparing 
the  amplitudes  of  the  data  with  the  synthetic  and  taking  the  average 
value.  In  general,  the  Inversion  yields  results  that  are  In  good 
agreement  with  the  literature.  He  will  briefly  discuss  the  Truckee 
event  as  an  example. 

The  Truckee  earthquake  was  a  strike-slip  event  at  10  km  depth 
which  produced  excellent  regional  records  but  very  few  teleselsmlc  body 
wave  records  as  typical  of  moderate  size  strike-slip  events.  The 

Truckee  earthquake  (a^  .  5.7)  ha$  been  studied  by  several  authors  (Tsai 
and  Akf,  1970,  and  Burdick,  1977)  making  It  a  good  test  case.  Tsai  and 
Akl  (1970),  from  first  motion  studies  and  modeling  of  the  surface  waves, 
determined  this  event  to  be  pure  strike-slip  on  a  fault  plane  striking 
N44°E  and  dipping  80°SE.  The  surface  wave  moment  was  determined  to  be 
25 

0.83  x  10  dyne-cm.  Figure  4.6  shows  the  location  of  the  epicenter, 

recording  stations  and  filtered  data  for  Truckee.  Also  shown  are  the 

synthetics  determined  from  the  Inversion  results.  Note  that  BOZ  and  TUC 

are  very  nearly  nodal.  The  Inversion  yields  a  mechanism  which  Is  very 

similar  to  Tsai  and  Akl's  (1970);  a  strike  of  N43°E,  a  dip  of  76°SE  and 

a  rake  of  -11°.  The  only  significant  difference  Is  the  slight  dip-slip 

component  In  our  solution,  which  Is  also  acceptable  on  the  basis  of  the 

first  motion  data.  The  moment  determined  from  the  waveforms  is  0.87 
2S 

xlO  dyne-cm,  which  Is  In  agreement  with  Tsai  and  Akl  (1970). 
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The  results  for  the  Orovllle  earthquake  are  particularly 
Interesting  In  that  the  observations  are  nearly  all  positive  as  opposed 
to  those  observed  te1ese1sm1ca11y(  see  Figure  3.5.  We  Mill  discuss  this 
event  later  after  Including  the  waveform  observations  In  the  12°  to  30° 
range. 


b)  Upper-mantle  responses 

At  distances  larger  than  about  12  the  mantle  "P"  wave  begins  to 
emerge.  This  feature  Is  displayed  In  Figure  4.9.  Near  1000  km  the 
beginning  portion  of  the  seismogram  Is  dominated  by  Pn,  that  Is,  a  P 
head  wave  traveling  along  the  top  of  the  mantle  as  a  head  wave.  Since 
this  earthquake  Is  nearly  •  pure  strike-slip  event  we  can  Interpret  the 
first  upswing  as  P^  followed  by  pPfl  and.sPn.  Beyond  about  12°  this 
combination  of  arrivals  becomes  much  sharper  since  they  are  no  longer 
head  waves.  Short-period  waves  from  explosions  fired  at  the  Nevada  Test 
Site  Increase  In  amplitude  by  nearly  two  orders  of  magnitude  near  this 
range,  for  example  see  Helmberger  (1973).  These  features  can  be 
generally  explained  by  the  tectonic  model  presented  In  Figure  4.10. 

This  smooth  model  Is  broken  Into  homogeneous  layers  suitable  for  the 
applications  of  GRT.  Essentially,  the  generalized  rays  returning  from 
the  Interfaces  In  the  mantle  are  added  to  the  rays  discussed  earlier 
(see  Figure  4.4).  Synthetics  appropriate  for  a  strike-slip  event  are 
presented  in  Figure  4.11.  The  recording  at  the  nearest  ranges  are 
essentially  the  same  as  displayed  In  Figure  4.5.  The  pattern  of  (P,  pP 
and  sP)  going  from  head  waves  to  diving  waves  Is  quite  apparent. 
Comparing  these  results  with  Figure  4.9  It  would  appear  that  the  lid 
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Fig  4.9  -  Profile  of  strike-slip  observations  from 
?nl  t0  Recordings  at  DUG,  ALQ,  and  LUB  ar«  from 
tne  22  December  1964  event  of  northern  Baja  Cali¬ 
fornia;  the  others  are  of  the  Truckee  event 
discussed  earlier. 


Fig  4.10  -  Velocity  depth  functions  for  T7 
(Burdick  and  Helmberger,  1978)  and  two  recent 
shear-wave  models  (Grand  and  Helmberger,  1982). 


Fig  4.11  -  (a)  Strike-slip  excitation  going 

from  Pni  domination  to  P  and  long-period  synthetics. 

(b)  Dip-slip  excitation  going  from  P^  domination 
to  P  and  long-period  synthetics.  v 

thickness  (the  high-velocity  layers  between  the  crust  and  the 

low-velocity  zone)  may  be  too  thin  In  the  T7  model.  Actually,  the 

crossover  distance  between  P.  to  P  shows  considerable  lateral  variation. 

n  - 

For  Instance,  short-period  Pn  In  shield  regions  are  observed  to  much 
larger  distances  than  In  other  regions.  However,  since  earthquakes 
rarely  occur  In  shield  regions,  we  will  restrict  our  discussion  to  model 
(T7)  which  appears  to  be  quite  useful  In  simulating  earthquake 
observations  In  active  regions,  see  for  example  Burdick  (1980). 


Nevertheless,  we  propose  omitting  records  at  ranges  12  to  14  as  being 
too  regional  or  structurally  dependent  to  be  useful  In  source  modeling 
in  most  situations. 

The  response  of  the  earth  models  at  ranges  beyond  15°  can  be 
obtained  by  simply  summing  primary  rays.  That  Is 
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The  term  "primary11  refers  to  rays  that  undergo  one  reflection  only.  The 
(a/slna)  factor  is  a  correction  for  spherical  spreading,  see  Gilbert  and 
Helmberger  (1972). 

Near  30°,  M(t)  approximates  a  step  with  amplitude  predicted  by  the 
curve  In  Figure  3.3.  Actually,  the  amplitudes  in  Figure  3.3  were  generated 
using  expression  (4.3)  assigning  the  JB  model.  Similar  results  are  obtained 
from  the  geometric  formula,  Bullen  (1965),  p.  126: 
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.  T  1$  the  travel  time 

where 

t  is  the  epicentral  distance  In  degrees 
10  Is  the  angle  of  Incidence 
vQ  Is  the  velocity  at  the  surface 
rQ  Is  the  radius  of  the  earth 

Because  the  ray  summation  has  omitted  multiples*  It  appears 
necessary  to  document  these  responses  by  comparison  with  other  methods. 
This  was  done  recently  by  Burdick  and  Orcutt  (1978)  with  a  comparison 
of  GRT  and  reflectivity  results  given  In  Figure  4.12.  These  same 
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Fig  4.12  -  (a)  The  slightly  smoothed  delta  function  responses 
from  a  reflectivity  calculation  (bottom)  are  compared  to  those 
from  a  generalized  ray  (top)  for  model  T7.  (b)  First  column 

contains  the  predicted  waveshape  due  to  the  source  only; 
second  column  Is  the  delta  function  response  of  the  model  T7; 
third  column  Is  the  convolution  of  the  first  two  columns. 

The  observed  waveforms  from  the  Borrego  Mountain  earthquake 
are  on  the  top  for  comparison. 
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responses  were  generated  by  "Full  Wave  Theory"  (see  Cormier  and  Choy, 
1981),  with  similar  waveforms. 

The  vertical  response  for  a  dislocation  source  can  be  generated 
from  M(t)  by  assuming  that  the  vertical  radiation  functions  Cj(p0)  etc., 
are  nearly  constant  over  the  (p's)  Involved,  accordingly 

<*•»  w  •  w*  j£  ffi  *  2  vm,i> 

0  L  j-i  J 

Similar  approximations  can  be  applied  to  generated  pP,  sP,  and  crustal 
multiples.  With  this  approximation,  the  procedure  becomes  Identical  to 
the  teleselsmlc  expression  (3.1  )  where  the  factor  (1/R)  Is  replaced 
with  a  convolution  of  M(t)  followed  by  a  derivative.  Synthetics 
appropriate  for  the  three  fundamental  faults  are  displayed  In  Figure 
4.13.  A  synthetic  seismogram  appropriate  for  ary  arbitrary  orientation 
can  then  be  formed  by  taking  a  linear  combination  of  these  three  as 
discussed  earlier. 

We  can  change  the  order  of  operations  and  generate  the  synthetic 
seismograms  as  If  the  M(t)  were  a  delta  function  as  under  the 
teleselsmlc  assumption.  The  first  column  of  Figure  4.12b  corresponds  to 
this  order  of  operations.  Convolving  with  the  H(t)  produces  the  proper 
synthetic  as  given  on  the  right  of  Figure  4.12b.  The  latter  order  of 
operations  Is  generally  used  for  studying  upper-mantle  structure  since 
M(t)  Is  Isolated  from  the  other  operations. 
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Fig  4.13  -  Upper-mantle  synthetics  without  and  with  Instrument  for  the  three 
fundamental  orientations  assuming  a  source  depth  of  8  km,  t*»l,  and  6ti  «  6t2  ■  6t»*l 
for  the  source  time  history. 


5.  Discussion 


In  this  review,  we  have  presented  a  procedure  for  generating 
synthetic  waveforms  for  bo<Jy  phases.  For  convenience,  we  divided  the 
eplcentral  distances  into  teleselsmlc,  regional,  and  local.  At 
teleselsmlc  distances,  we  showed  that  long  period  boty  phases  from 
shallow  earthquakes  are  coherent  at  neighboring  stations  and  that  the 
observed  waveform  could  be  decomposed  In  a  manner  that  allows 
determination  of  faulting  parameters.  It  appears  that  mar\y  events  In 


the  Magnitude  class  5.5  to  6.5  can  be  easily  Modeled  at  long  periods  by 
assuming  point  sources.  This  result  should  not  be  surprising  since  the 
signal  wavelengths  are  long  compared  to  the  fault  dimensions.  On  the 
other  hand,  the  corresponding  comparison  between  short  period  records  at 
neighboring  stations  Is  much  less  coherent  and  complicated,  see  for 
example  figure  5.1.  Some  success  at  explaining  these  records  has  been 
presented  recently  by  Ebel  and  Helmberger  (1982)  by  Introducing  complex 
faulting. 

To  model  shorter  periods  requires  that  we  generalize  our  point 
source  procedures  to  treat  finite  faults  as  distributed  point  sources, 
see  for  example  Langston  (1978)  and  Heaton  (1982).  Following  the  latter 
technique  we  simply  simulate  rupturing  on  a  fault  plane  by  summing  point 
sources  with  the  proper  timing,  see  figure  5.2.  The  procedure  Is  made 
somewhat  complicated  by  the  presence  of  the  free  surface  which  requires 
summing  pP,  etc.  Thus,  teleselsmlc  signals  can  show  considerable 
variation  with  azimuth.  Note  that  the  long  period  records  prove  quite 
Insensitive  to  the  rupture  properties  In  this  simple  case  compared  with 
the  short  pelrfods.  Conceptions) ly,  we  should  be  able  to  determine  the 
time  structure  of  the  rupture  process  by  simply  matching  short  period 
waveform  patterns  and.  Indeed,  this  subject  Is  being  pursued  vigorously. 
Unfortunately,  at  these  periods  the  many  Idealistic  assumptions 
Involving  the  homogeneity  of  the  mantle,  simple  flat  layering  near  the 
source  and  receivers,  planar  fault  planes,  to  name  a  few,  begin  to  break 
down.  Thus,  the  pattern  recognition  Is  not  so  obvious. 

A  simpler  approach  is  to  suppose  that  the  rupture  process  can  be 
simulated  with  some  small  number  of  point  sources  where  each  source  has 
Its  particular  time  history.  By  modeling  the  long  and  short  periods 
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Fig  5.2  -  Comparison  of  the  short-period  and-  long-period  synthetics 
at  the  teleseismic  stations  for  various  fault  models.  The  numbers 
are  the  amplitudes  In  micrometers  calculated  for  the  seismic  moment 
of  .5  x  1024  dyne-cm.  The  top  row  corresponds  to  a  particular 
distribution  of  faulting  appropriate  for  the  Coyote  earthquake 
where  the  dislocation  Is  largest  at  the  center  and  tapers  down  at 
the  edges.  The  Ul,  U2,  U3,  and  U4  models  show  uniform  slip 
histories  but  with  different  rupture  speeds  (After  Liu  and 
Helmberger,  1982). 


simultaneously  with  these  few  parameters,  we  can  look  for  sub-structure. 
Such  a  procedure  was  tried  recently  by  Wallace  el  el  (1981b)  -on  the  Santa 
Barbara  earthquake  of  1978,  see  figure  5.3. These  results  Indicate  that 
on  a  smaller  scale  most  earthquakes  are  Indeed  complex  as  suggested  by 
Kanamorl  and  Stewart  (1978).  The  sharp  jumps  Indicated  in  the  time 
history  of  figure  5.3  are  obviously  Important  with  respect  to 


Fig  5.3  •  Modelling  experience  obtained  from  studying  the  near-  and  far- 
fleld  records  of  the  1978  Santa  Barbara  earthquake,  a)  displays  the 
far-fleld  time  history  for  an  effective  point  source;  b)  displays  the  fit 
of  the  long-period  synthetics  and  data;  c)  displays  the  short  period  fit. 
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Fig  5.5  -  Comparison  of  synthetics  with  waveform  data  at  all  ranges 
with  preferred  model  (e  ■  215  ,  X  ■  -65,  and  6  ■  48°).  Inversion 
results:  with  the  5  P*4  records  exclusively  are:  (e  »  195°, 

X  *  -72°.  end  6  ■  46°TT  10  upper-mantle  ranges  exclusively  are: 

(e  »  197°,  X  ■  -63°,  and  5  »  58°),  8  teleselsmlc  waveforms  exclusively 


ere:  (e  •  221 


and  S  ■  44  )  (After  Yao  et  el.,  1982). 


understanding  the  faulting  mechanism,  but  their  existence  Is 
particularly  significant  In  estimating  earthquake  hazards.  That  Is,  the 
high  frequency  strong  motions  appear  to  be  more  strongly  controlled  by 
the  jumps  in  the  time  history  than  by  the  overall  duration  (see  Wallace 
et  a1_.,  1981b;  Ebel  and  Helmberger,  1982). 

At  ranges  less  than  about  30°,  we  found  that  the  bocfy  wave 
arrivals  become  multi-valued  and  some  care  needs  to  be  taken  to  remove 
propagatlonal  features  from  source  phenomenon!.  Some  progress  In 
achieving  this  separation  can  be  seen  In  figures  5.4  and  5.5  where  we 
show  that  consistent  estimates  of  the  source  orientation  can  be  obtained 
at  nearly  aqy  range.  Note  that  there  are  no  teleselsmlc  waveforms  for 
the  Truckee  event  compared  with  the  Orovllle  event,  even  though  they  are 
roughly  the  same  strength.  The  reason  Is  due  to  the  difference  In 
vertical  radiation  pattern  as  discussed  earlier  (see  figure  3.2).  How 
successful  this  analysis  will  be  for  other  regions  Is  not  known  but 
Initial  results  look  promising. 
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Appendix  A.  Generalized  Ray  Theory 


The  method  of  generalized  ray  theory  allows  one  to  construct 
the  solution  of  a  seismic  boundary  value  problem  by  Inspection.  The 
method  was  first  Introduced  to  the  seismic  community  by  Spencer  (1960). 
The  Idea  can  be  easily  understood  by  working  through  a  simple  example, 
namely  Me  sandwiched  fluid  layer  problem  (see  figure  Al).  The  solution 
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with  the  simple  boundary  conditions  on  continuity  of  vertical  stress  and 
displacement  across  the  Interfaces,  that  Is 


These  four  equations  are  used  to  solve  for  the  four  unknowns  A,  B,  C, 
and  D.  Substituting  the  fj  Into  the  boundary  conditions,  we  obtain 
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Each  value  of  (n)  corresponds  to  a  particular  ray,  for  example 
with  n  •  1  we  obtain 
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where  the  Identity  (1-R^y  )  •  we*  used. 
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and,  similarity,  for  n  ■  2  we  obtain 
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Thus,  the  solution  Is  built-up  by  adding  up  the  various  rays  which  can 
be  determined  by  Inspection. 

He  can  check  our  solution  and  resulting  code  by  generating  the 
step  response  of  a  sanArtched  layer  In  the  limit  of  small  thickness  and 
compare  the  results  with  the  solution  with  no  layer,  see  figure  A2. 

Next,  we  examine  the  case  where  we  go  beyond  critical  angle  so  that  the 
receiver  Is  In  the  geometric  shadow  and  Investigate  tunneling.  In  this 
case,  the  layer  Is  too  thin  to  support  a  head  wave  and,  thus,  the  strong 
Interference  effects  Introduced  by  the  multiples.  The  solution  can  be 
written 
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expanding  D  In  a  series  yields  the  ray  expansion  and  the  multiples 
displayed  In  figure  A2.  When  the  layer  Is  too  thin  to  support  head 
waves,  we  can  obtain  a  useful  approximation  by  noting  that 

Limit  D  -  T-  -  -  1 
Th  *►  0  X 

Thus,  tht  only  effect  of  the  layer  on  the  solution  Is  In  the  contour 
t  •  pr  4^0*  4  H)  +  Th  tij 

Essentially,  the  contour  moves  off  of  the  real  p-axls  at  p  »  l/ag 
and,  thus,  no  longer  will  there  be  a  square-root  singularity  at  the 
direct  arrival  time  but  slightly  smoothed  depending  on  the  thickness. 
That  Is,  the  height  of  the  steps  given  In  figure  A2  corresponds  to  the 
same  amplitude  of  the  steps  with  no  layer.  However,  the  rise- time  of 
these  steps  are  no  longer  sharp.  Thus,  when  these  step  responses  are 
used  to  generate  synthetics,  they  filter  out  the  high-frequencies  or 
short  periods.  The  effect  Is  nearly  Identical  to  the  Q-operators 
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Fig  A2  -  Numerical  step  response  of  tunneling. 
The  source  depth  Is  3  km  with  a  20  km  source- 
receiver  separation.  Seismic  parameters  are 
Pi«  2.3  gm/an3,  Bi  *  2.02  km/sec,  and  p2«  2.4 
gm/cm3,  e2«  2.31  km/sec,  with  a  barrier 
thickness  of  10  meters.  Relative  peak 
amplitudes  are  given  for  each  case,  n  Indicates 
the  number  of  rays  Included  In  the  sum. 


discussed  earlier. 

The  technique  of  constructing  a  solution  by  simply  adding  up 
rays  that  have  traversed  the  various  layers  with  some  specified  mode  of 
propagation  Is  called  the  method  of  Generalized  Reflection  and 
Transmission  Coefficients."  The  coefficients  for  the  solld/solld 
Interface  are  somewhat  more  complicated  because  of  the  possible  mode 
changes  from  P  to  SY(  etc.,  but  the  procedure  Is  basically  the  same. 

The  generalized  reflection  and  transmission  coefficients  used  In  this 
paper  are  given  by  Helmberger  (1968). 

Another  concept  commonly  used  In  GRT  Is  that  of  receiver 
functions.  Note  that  we  have  been  working  with  potentials  throughout 
these  notes.  Essentially,  we  start  with  a  source  potential  and 
propagate  It  to  some  location  by  applying  the  proper  product  of  T  and  R 
coefficients.  At  the  receiver  we  need  to  convert  this  potential  Into  a 
displacement  or  some  measurable  quantity*  In  a  whole  space,  the 
displacements  can  be  obtained  by  simply  performing  the  derivatives 
defined  by  equation  (2.24).  On  a  free  surface  the  conversion  from 
potential  to  displacement  becomes  more  complicated,  see  figure  A3.  The 
proper  way  to  satisfy  the  boundary  conditions  Is  to  take  the  limit  as  Z 
goes  to  zero  and  sum  the  three  arrivals  since  the  travel  paths  coincide 
and  label  the  results  "receiver  functions."  The  results  for  Incoming 
P-waves  art: 

*p2  ■  2n0(ns2  -  p2)/e2»(p) 

RpR  •  “4V>gP/02i>(p) 

whara 
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Fig  A3  -  Schematic 
showing  the  Interaction 
of  the  three  arrivals 
that  form  the  receiver 
function. 


The  second  subscript  Indicates  component,  namely  Z  for  vertical  and  R 
for  radial.  The  results  for  Incoming  SV-waves  are: 


RSZ  "  HV82d(p) 

*sr  "  2T,e(V*“  ®2  D(p) 


Table  1 
Crustal  Model 

SVEL  Density  Layer  Thickness 

3.5  2.7  32.0 

4.5  3.4 


PYEl 

6.2 

8.2 


